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Abstract 

Computational tools for normal mode analysis, which are widely used in physics 
and materials science problems, are designed here in a single package called NMscatt 
(Normal Modes & scattering) that allows arbitrarily large systems to be handled. 
The package allows inelastic neutron and X-ray scattering observables to be calcu- 
lated, allowing comparison with experimental data produced at large scale facilities. 
Various simplification schemes are presented for analysing displacement vectors, 
which are otherwise too complicated to understand in very large systems. 
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1 Introduction 

At large scale facilities for neutron and X-ray scattering, large quantities of ex- 
perimental data are produced. For complex, nanoscale systems, understanding 
this data requires computer models. In the case of inelastic scattering, molec- 
ular dynamics (MD) simulations [1] are widely used to equilibrate structures 
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and explore dynamics as a function of temperature and other experimen- 
tal parameters. However MD only gives a partial description of vibrational 
modes through the partial density of states and when knowledge about spe- 
cific vibrational modes is required, normal mode analysis (NMA) [2] has to 
be performed. For physics and materials science problems, NMA gives a de- 
scription of the lattice dynamics via the dispersion (k-vector dependence) of 
the mode frequencies [3]. For small systems (< 200 atoms in the simulation 
box) very accurate results can be obtained using density functional theory 
(DFT) methods to determine interatomic force constants [4,5]. By combining 
DFT and software that constructs and diagonalises the dynamical matrix and 
calculates the experimental observables, experimentalists now have sophisti- 
cated tools to analyse their data. The PHONON [6] package is one of the best 
examples. 

Phonon codes are traditionally limited to small systems for a number of rea- 
sons. For example, for small unit cells, the reciprocal lattice is big and stronger 
effects of dispersion are expected. If DFT methods are used to determine 
force constants then these methods are themselves restricted to a few hun- 
dred atoms. However the development of nanoscale structures of (partially) 
crystalline materials stimulates a need for phonon codes to be extended to 
much larger systems. Parameterised force fields [7] can be used to determine 
the inter-atomic force constants and, as will be seen in the example presented 
here, strong dispersion effects are observed inspite of the small reciprocal lat- 
tice. 

In biomolccular systems, the need for NMA has long been recognised and codes 
like CHARMM [8] allow the gamma point normal modes to be calculated for 
moderately big systems. In addition, the neutron scattering quantities can be 
directly calculated from the simulations using the time-correlation function 
formalism [9], as implemented in the nMOLDYN program [10]. A combina- 
tion of neuton scattering experiments and atomic detail computer simulations 
has proven to be a powerful technique for studying internal molecular vi- 
brations [11,12,13]. In this approach one can validate the applied numerical 
models, i.e. force field parametrizations, depending on the agreement between 
the experimental and calculated spectra. 

In this paper we present a software package that extends the functionality 
of codes like PHONON [6] and Climax [14] to arbitrarily large systems and 
extends the gamma point only analysis already available for larger systems 
to include k-vector dependence. The software reads a Hessian matrix of force 
constants, constructs and diagonalises the dynamical matrix for any k-vector 
and calculates neutron and X-ray scattering observables. The computational 
bottleneck remains the diagonalisation of correspondingly large dynamical ma- 
trices and we comment on approximations that have to be used when the Bril- 
louin zone cannot be sampled at a large number of points. In large systems. 
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atomistic detail in the displacement vectors can be difficult to interpret due 
to the large number of degrees of freedom and we present two methods for 
simplifying this information. The first entails summing displacement vectors 
over atoms in user-defined beads, while the second involves a reduction of the 
degrees of freedom in the dynamical matrix by summing over force constants, 
which has the advantage of reducing the number of modes to be examined. 



2 Theoretical background 



The standard approach, also called a direct method [15], to the lattice vibra- 
tion problem of crystals is based on the explicit knowledge of the interaction 
between all atom-pairs in the system. Subsequently, one deduces the corre- 
sponding force constants, and constructs and diagonalizes the dynamical ma- 
trix for any k- vector in order to obtain the frequencies of the normal modes. A 
reasonable atomic detail description of interactions within large biomolecular 
systems are provided using empirical force fields. 

In the following we will briefly summarize the aspects of the classical theory 
of lattice vibrations [16] and proceed to the description of the explicit phonon 
calculations. 

The individual atomic positions in the crystal can be assigned as 

Rnix{t) ^Rn + r^ + Unij.{t), (l) 

where Rn is unit cell lattice vector and Un^{t) is displacement of atom // from its 
equilibrium positions r^. Within the harmonic approximation we concentrate 
on expansion of the small differences of potential energy V due to the small 
changes in atom positions: 

1 

V(u) pa Vo + X] ^ ^"Z*" +9 Un^J^n^ia,mu(iUmuj3 + ■ ■ ■ , (2) 



where the second derivative defines the force constant between the atoms /j, 
and u: 

'Dnu,a,mwP = ^ ^ , («, P^X,y, z) (3) 



As each unit cell is identical to every other unit cell in the crystal, the dis- 
placement pattern of a normal mode has to be identical to that in any other 
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cell to within a phase difference k{Rn — Rm)- The representation of the atom 
displacement is chosen to be a plain wave ansatz of the form: 



u 



nfj, 



{k,t) 



ei,k^Mi[kRn-uJj:t\), 



(4) 



where e^^ is the polarization vector and is the mass of the atom /x. We 
omit writing Cartesian component subscripts. Solving the equation of motion 
with ansatz (4) is equivalent to the eigen-value problem 



where 



^M'^ (^) = E / ^ ^nn,mu eX.p[ik{Rm - R„)] (6) 



is the so called dynamical matrix. 

The form of the dynamical matrix (6) requires the atom pairs for which one 
atom belongs to a different unit cell, i.e. (m ^ n), to be identified. These 
terms contribute to the so called Bloch-factor cxp[ik{R„i — Rn)] and make the 
dynamical matrix complex. But in case of applying periodic boundary condi- 
tions (PEG) as implemented in computer simulation programs, the potential 
energy of a crystal is given as an explicit function of only the atom positions in 
the primary unit cell. As a consequence, we obtain the second derivative ma- 
trix H in which the contributions from the inter-cell atomic pairs are mapped 
and added to the corresponding image atom pairs in the primary unit cell: 



One can directly obtain 'Dn^ia^mup by increasing the size of the unit cell by 
one or more layers of periodically arranged image cells and calculate force 
constants in the extended supercell. However, this approach is unfavourable 
when dealing with very large systems. 

A similar approach is to decompose the potential energy V^^'-^ in equation (7) 
into individual contributions from the image cells, V^^^ — Vq -\- X^m Kn and 
evaluate the second derivative matrix for each term upon the same minimized 
structure. 

The situation is simpler if the interaction is truncated at some cutoff distance 
Tcut so that the "minimum image convention" (MIC) is obeyed. The MIC 
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states that each atom interacts at most with one image of every other atom in 
the system (which is repeated to fully enclose the primary unit cell with the 
periodic boundary conditions). This has the effect of limiting the interaction 
cutoff, for example, to no more than half the length of the minimum side 
when simulating the orthorhombic cell, Vcut < min{a/2,b/2, c/2}. It should 
be noted that the size of nanoscale crystals usually far exceeds the spatial range 
of forces between atoms (~ <12 A) allowing physically reasonable cutoff radii 
to be introduced. 
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Fig. 1. A 2D periodic system representing the effect of the minimum image con- 
vention while constructing dynamical matrix. If a distance between atom v (hollow 
blue circle) and fi (full black circle) is larger than the cut-off distance (red circle 
around the atom fi), and if the force constant H^u turns out to be nonzero, it means 
that the value of H^u has to be related to the atom i' in the image cell m (full blue 
circle), with the position vector fv + Ix (Ix = Rm — Rn), instead of to the atom 
fv in the primary cell n. The corresponding lattice translation vector Ix is thus 
determined for the element H^,y. 

According to the MIC we can conclude that there is always only one translation 
t per atom pair (i/, fi) giving rise to the minimum distance |(rv + — r^\- 



— { — 1, 0, l}a=i,2,3 • t — Sill + S2I2 + 53/3; 



so that: |(tv + t) 



mm. 



where la are lattice translation vectors, c.f. Fig. 1. In the case of a nonvanishing 
H^^i, we get the following expression for a dynamical matrix element 

- 1 I H^,i^> ii:\r^-r^\<rcut 

T^A^) = r—— { (9) 
'M^M^ H^,^exp(iH), if : - > r^ut 



According to equation (5) the diagonalization of matrix 'D{k) yields the phonon 
frequencies tu^^- and corresponding polarization vectors e^^^. for a given phonon 

wave vector k. A complete solution leads to the phonon dispersion relations. 
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The sTibscript j denotes a branch in the phonon dispersion. In a crystal of N 
atoms, there are 3A^ branches. 



Inelastic scattering 

The dynamic structure factor S{q, uj) contains information about the structure 
and dynamics of the sample. It can be split into a coherent part arising from 
the cross-correlations of atomic motions and an incoherent part describing 
self-correlations of single atom motions. According to the standard theory 
[17,18], which is based on the harmonic approximation, we obtain the following 
expressions for the coherent and incoherent dynamical structure factors: 



G kj ^"^h 



'-'^^■''^^^^■exp(-W^^(^)+igr-) 




X 



kj) 



(10) 



and 



'5(g,u;),nc = E'^rE ^T^k"- v/(^^fe,-) + 1) X 

xexp[-2W^(g-)]5(cu-^^.), (11) 

where q is the scattering vector, 0"^^ is the corresponding atomic scattering 
length, [n[iS) -|- 1) refers to the phonon creation process (absorption spectrum) 
and n[uS) is the mean number of phonons of frequency uj at temperature T 
according to the Bosc- Einstein statistics 

exp(/.^/A:^T) 



The factor exp[— 2FF^(g)] is called the Debye- Waller factor: 
-» g/ \ -» 

W^m(^) = 7, > ^ocM = («/xa'"^^)T, (13) 



where B(/x) is a 3 x 3 symmetric tensor representing the thermodynamic mean 
square displacement of an atom which can be expressed by the partial 
atomic phonon density of states ga^^^iLo): 



^"^^^^ = 2M; I ^3a,,M coth J . (14) 
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The partial atomic density of states is a weighted distribution of normal modes 

2 n3N 



where n is the number of sampling ^-points in the first Brillouin zone. 

The evaluation of the Debye- Waller B factors using equation (14) requires 
extra attention due to the "zero-phonon" term resulting from the singularity 
caused by phonon (acoustic) branches where u){q,j) = 0. The contributions of 
the 3 acoustic modes is treated separately by using the Debye approximation 
for the density of states, i.e. Qad^) oc cj^, and normalization jg^™"^ gac{uj)duj = 
3, where Umax is the maximum frequency up to which the acoustic dispersion 
curve is linear 



jjj max . 

i \-sr f * f I - 



j=l ^'^fi"] max Q 
j=l ^jx^j max 




I 



keT y f X , 1 

-dx + 



fl^jmax) J exp(x) - 1 4 



(16) 



A new variable Xjmax — T^jmax/ksT was introduced in the last equation. 



3 Analysing the displacement vectors 



For systems containing thousands of atoms (N) the displacement vectors ob- 
tained by diagonalising the dynamical matrix can be difficult to understand, 
especially for low frequency modes which involve the displacement of many 
(or all) atoms. One simple solution to this problem is to sum over the dis- 
placements of atoms within beads, which represent logical coarse grains of the 
system, for example base molecules in the DNA example below. This treatment 
allows different bead definitions to be applied to the calculated displacement 
vectors but has the disadvantage of not reducing the number of displacement 
vectors from 3N. 

A related approach is to reduce the atomic level Hessian matrix to lower 
dimension by mapping the inter-atomic force constants on to inter-bead force 
constants [19]. For N' beads, the resulting dynamical matrices have dimension 
3N' and therefore result in 3N' displacement vectors for any k-vector. We 
note that any reduction in the dimensionality of the system causes a loss in 
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Fig. 2. Coherent dynamical structure factor of B-DNA. Scattering vector qis varied 
along the helical axis (0,0,1). Shown are S{q,u;), \q\ = 0.019 . . . 0.5811"^ 

information, which is the rotational degrees of freedom of the beads (rigid 
bodies). 



4 Example 

To verify the implemented formalism we have simulated a B-form DNA molecule 
(right-handed, 10 base-pairs per turn, pitch 33. 6A) using CHARMM [8]. The 
full crystal environment was generated using periodic boundary conditions for 
an orthorhombic unit cell containing one helix of DNA. The dimension of the 
unit cell is a = 32. 2A, b = 31. 8A and c = 33. SA, with c parallel to the helical 
axis. The starting configuration was obtained by minimization of the potential 
energy of the crystal structure obtained as the time average over a Ins MD 
simulation at lOOK. The Hessian matrix of force constants was generated by 
displacing each atom in turn from equilibrium and calculating the forces in- 
duced on all other atoms. Diagonalisation of the resulting dynamical matrices 
was performed using the the routine zcheev from the LAPACK library [20]. 

The typical coherent spectrum of DNA, Figure 2, obtained from equation 
10, shows a well-defined Brillouin peak at small uj which moves along the 
frequency dimension upon varying momentum transfer q. Fitting the spectral 
profile with a Gaussian as a function of wave- vector gives the dispersion curve 
shown in Figure 3, which compares well with the recent experimental results 
[21]. Figure 3 also shows the result of the equivalent analysis of a 300K MD 
simulation on the same model of B-DNA using nMoldyn. 
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Fig. 3. Dispersion curves obtained by Gaussian fitting of the calculated Brillouin 
peak for inelastic X-ray scattering using NMscatt (solid curve + circles), and using 
nMoldyn (dashed curve + squares). 



In order to gain insight into the nature of low frequency dynamics we can 
analyse the displacement vectors at an atomic level (see Figure 4(a) for an 
acoustic mode). By summing over the displacement vectors in terms of beads, 
where base molecules, sugar molecules and phosphate groups are treated as 
single units, a simplified picture of the normal modes is obtained. Figure 4(b) 
shows a high frequency mode which has a pronounced contribution from the 
phosphate groups. 



5 Conclusion 



The new, user-friendly computational package NMscatt presented here enable 
an efficient atomic detail analysis of different types of inelastic scattering ap- 
plied to arbitrarily large nanoscale systems. The ability to perform molecular 
dynamics and phonon calculations on large nano and bio-molecular materi- 
als means that one can efficiently pursue the investigation of some poorly 
understood structural and dynamical features of these systems. 
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(a) B-DNA: acoustic mode (b) B-DNA: beads 

Fig. 4. Atomic a) and h) beads (base molecules, sugar molecules and phosphate 
groups) displacements of the selected mode. 
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A Program package and data structure 

There are four main modules in the NMscatt program package phonon, coh, 
incoh and bead, and the overall NMscatt structure is given in Figure A.l. Below 
are described the corresponding modules. 

• phonon: Providing the full Hessian matrix for a given energy-minimized 
atomic structure within the specified crystallographic unit cell this module 
constructs dynamical matrix and calculates its eigenvalues and eigenvec- 
tors at given wave vector k. At input this module requires to specify the 
Bravais lattice vectors that were previously used in the molecular mechan- 
ics/ dynamics simulation package to satisfy the periodic boundary conditions 
while generating the minimized structure. The cut-off radius must be given 
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at which the long range interactions arc truncated while calculating the Hes- 
sian matrix in the simulation. At output two separate files eig_val_n and 
eig_vec_n are written containing eigenvalues and eigenvectors, respectively, 
where n assigns a consecutive number of the sampling fc-point in the first 
Brillouin zone. These files serve as an imput for other modules of NMscatt. 

— * 

The choice n = is assumed to be reserved for the F-point, k = (0, 0, 0). The 
lowest n should sample the vicinity of the F-point. The /c-points are to be 
specified in the fractional coordinates with respect to the reciprocal lattice 
vectors. The elements of the upper triangle of the hessian matrix should be 
provided in the binary (default) hessian_uf or alternatively in the ASCII 
file hessian_f and the atomic coordinates in the file coord writen in the 
CHARMM coordinate format. 

• incoh: This module allows to calculate atomic Debye- Waller factors repre- 
senting the mean square displacements and dynamical structure factor of 
incoherent one-phonon neutron scattering on monocrystals and from orien- 
tationally averaged powder. In the latter case we need to specify an abso- 
lute value of momentum transfer q, range of the /c-point index n: — rimax 
for picking up the corresponding A;-point eigenvalues- and eigenvectors-files 
generated by phonon, number of random orientations of q vector to provide 
spherical averaging and absolute temperature. Also, one has to define the 
/c-point index Ua, for which the /c-point range — corresponds to the 
linear regime of the acoustic-mode-dispersion curves. This is needed for the 
proper derivation of Debye- Waller factors using Debye approximation. As 
a result, the function S{q,u) is given in the file s_q_w. Optionally, one can 
also obtain density of states (DOS) in this module. 

• coh: This module calculates dynamical structure factor of coherent one- 
phonon neutron or X-ray scattering on monocrystals. In addition to the 
input data required by incoh (except for the spherical averaging), we need 
to define for module coh also the type of scattering neutron/X-ray, the 
number of higher Brillouin zones included for sampling momentum transfer 
vector, and the range of the A:-point index ni — n2 evaluated by module 
phonon and assigning the fc-points in the first Brillouin zone, which lie along 
the selected direction of the momentum transfer vector q. As a result, the 
function S{q, ur) is given in the file s_qw_coh. 

• bead: This module is used to enable visualization of selected vibrational 
modes obtained by running phonon for the F-point, such that atomic dis- 
placement vectors are projected on to the beads, which arc defined as the 
center of mass of larger atomic groups of the system, for example residues. 
The output files are readable by the program xmeLkeinol[22] which enables 
direct visualisation of the mode. 
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CRYSTAL STRUCTURE 
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Fig. A.l. Flow-chart of the NMscatt program package. 

Compiling 

To compile the program package Makefile is provided. It is important to note 
that 64-bits processors are prerequisite for applying the NMscatt to analy- 
sis of the large systems (containing more than 2000 atoms). One needs to 
install the LAPACK library on the computer beforehand, prior to NMscatt. 
The fortran compiling switches g77 -mcmodel=medium -f unroll-all-loops 
-fno-f 2c -03 are recomended when installed on 64-bits processors running 
Linux. 



Benchmark results 



Bencmark results were obtained on AMD Athlon 64 X2 Dual Core Proces- 
sor 2.2GHz running Linux for B-form DNA simulated with the CHARMM 
program. 
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